Dissipative hydrodynamics for viscous relativistic fluids 
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Explicit equations are given for describing the space-time evolution of non-ideal (viscous) rel- 
ativistic fluids undergoing boost-invariant longitudinal and arbitrary transverse expansion. The 
equations are derived from the second-order Israel-Stewart approach which ensures causal evolution. 
Both azimuthally symmetric (l-l-l)-dimensional and non-symmetric (2-|-l)-dimensional transverse 
expansion are discussed. The latter provides the formal basis for the hydrodynamic computation of 
elliptic flow in relativistic heavy-ion collisions including dissipative effects. 
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I. INTRODUCTION 

Ideal fluid dynamics has been used successfully to pre- 
dict the collective flow patterns in Au-|-Au collisions at 
the Relativistic Heavy-Ion Collider RHIC (for a review 
see 0). The ideal fluid description works well in al- 
most central Au+Au collisions near midrapidity at top 
RHIC energy, but gradually breaks down in more periph- 
eral collisions, at forward rapidity, or at lower collision 
energies 0, indicating the onset of dissipative effects. 
To describe such deviations from ideal fluid dynamics 
quantitatively, and to use the experimental data for the 
extraction of values or phenomenological limits for the 
transport coefficients of the hot and dense matter cre- 
ated during the collision, requires the numerical imple- 
mentation of dissipative relativistic fluid dynamics. Al- 
though a formulation of such a theory which avoids the 
longstanding problems of acausal signal propagation and 
other instabilities associated with the original relativis- 
tic fluid equations given by Eckart jSj] and Landau and 
Lifshitz has been known for almost 30 years sig- 
niflcant progress towards its numerical implementation 
has only been made very recently 0, S S ■ this 
point, we are only at the very beginning of a program 
that will eventually apply viscous relativistic fluid dy- 
namics to heavy-ion collision data. Existing numerical 
implementations are (l-l-l)-diniensional and can only de- 
scribe cylindrically symmetric transverse expansion with 
boost-invariant longitudinal dynamics 0, |^ . As we will 
show here, even the (l-l-l)-dimensional case still presents 
some open formal issues which we address in the present 
paper. The numerical codes are still in the process of 
being tested and will not be discussed here. 

The paper is organized as follows: In Section ^1 we 
shortly review relativistic ideal fluid dynamics and the 



conditions for its applicability. While most of this is stan- 
dard textbook material, it helps to establish notation and 
to better appreciate the differences in the non-ideal case. 
In Section UTTI we discuss the non-ideal fluid decomposi- 
tion and introduce the dissipative flows (bulk and shear 
viscous pressure, heat conduction) and how they manifest 
themselves in the baryon current and energy-momentum 
tensor in the Eckart and Landau frames. Section llVl deals 
with the derivation of equations to determine the evolu- 
tion of these dissipative flows. We follow the treatment of 
Israel and Stewart Q and discuss both the (acausal) flrst- 
order and (causal) second-order theories (this nomencla- 
ture will be explained in Sec. llVjl . While most of the 
material up to this point can already be found elsewhere 
(see also the beautiful lecture notes by Rischke in 
jlOl|). it is needed here for a selfcontained presentation 
and for a critical discussion of some systematic expansion 
issues which we point out in Sec. IIV Bl and which are of 
practical relevance. Section ]V\ contains the main results 
of this paper (with many technical details deferred to the 
Appendix) , namely complete sets of causal equations of 
motion for the dissipative transverse hydrodynamic ex- 
pansion of systems undergoing boost-invariant longitu- 
dinal flow. The discussion of the azimuthally symmet- 
ric (l-l-l)-dimensional case in Sec. IV Al improves on the 
presentation given in the recent work by Muronga and 
Rischke 8], while the equations for the non-symmetric 
(2+l)-dimcnsional case in Sec. IV Bl are original and have, 
to our knowledge, not been presented before. The con- 
cluding Section lVll summarizes our results and gives some 
further discussion. 



II. IDEAL FLUID DYNAMICS 
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Before explaining the structure of the equations for 
causal dissipative relativistic fluid dynamics, let us 
quickly review the case of ideal fluid dynamics. Any fluid 
dynamical approach starts from the conservation laws for 
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the conserved charges and for energy-momentum, 



d^Tt"" = 0. 



0, z = l,. 



1 



(2.1) 
(2.2) 



For simpHcity we will restrict ourselves to fc = 1 (say, 
Nf^ = net baryon number current) and drop the index 
i in (|2.1|l . It must also ensure the second law of thermo- 
dynamics 



9^5^ > 0> 



(2.3) 



where S'^ is the entropy current. Ideal fluid dynamics fol- 
lows from these equations under the assumption of local 
thermal equilibrium, i.e. if the microscopic collision time 
scale is very much shorter than any macroscopic evo- 
lution time scale such that the underlying phase-space 
distribution f{x,p) relaxes essentially instantaneously to 
a local equilibrium form (upper signs for fermions, lower 
signs for bosons) 



/oq(a;,p) 



1 



^\p-u{x)+iJ.{x)\/T{x) J^l' 



(2.4) 



where u^{x) is the local fluid velocity at point x, fJ-{x) 
is the local chemical potential associated with the con- 
served charge N (it enters with opposite sign in the dis- 
tribution / for antiparticles), and T{x) is the local tem- 
perature. Plugging this into the kinetic theory definitions 



N^ix) 
T^'^'ix) 
S^{x) 



^p^p'^[f(x,p) + fix,p)l 



,f{x,p) ln/(a;,p) 



±(lT/(x,p)) \n{lTf{x,p)) + {f 



(2.5) 
(2.6) 
(2.7) 

/)", 



leads to the ideal fluid decompositions 



m^^nu^, (2.8) 
T^^^ = eu^w'' -pA'"' (with A'"'=g'"'-u''u''), (2.9) 
^e^q^su", (2.10) 

where the local net charge density n, energy density e, 
pressure p and entropy density s are given by the stan- 
dard integrals over the thermal equilibrium distribution 
function in the local fluid rest frame and are related by 
the fundamental thermodynamic relation 



T s = p — iin 



(2.11) 



Inserting Eqs. ^T^-^TT^ into Eqs. (EHJ and yields 
the relativistic ideal fluid equations shown in Eqs. (|2.12|) - 
(|2.14() below. Using Eq. (|2.11() together with the Gibbs- 
Duhem relation dp — s dT + ndfj,, it is easy to prove that, 
in the absence of shock discontinuities, these equations 
also conserve entropy, i.e. d^S'^ = 0. 



Note that the validity of the decompositions (|2.5|) - 
(|2.7|) only requires local momentum isotropy (i.e. that 
in the local fluid rest frame the phase-space distribu- 
tion reduces to a function of energy E only, f(x,p) = 
f{p-u{x);T{x), iJ,{x)^), but not that the distribution 
function has the specific exponential form H2.4|) which 
maximizes entropy. This may have relevance in situations 
where the time scale for local momentum isotropization 
is much shorter than for thermalization 0, 0, 0| (i.e. 
it is much easier to change the direction of the particles' 
momenta than their energies), with the macroscopic hy- 
drodynamic time scale in between. In this case the lo- 
cal microscopic states would not maximize entropy, and 
the relation (|2.11|) would not hold between the quantities 
e, p, n, and s defined through eqs. (|2.5|) - (|2.1U|) . Still, they 
would follow ideal fluid dynamical evolution since en- 
tropy production by microscopic kinetic energy-shifting 
processes would happen only on time scales which are 
large compared to the macroscopic evolution time scales. 
[Note that in the absence of such a clear separation of 
time scales entropy production can not be neglected dur- 
ing the macroscopic evolution, and ideal fluid dynamics 
must be replaced by dissipative fluid dynamics.] 

The ideal fluid equations read (with 9 = d-u denoting 
the local expansion rate) 



n 

e : 



-n9, 
ie+p)( 
Vp 
e+p' 



(2.12) 
(2.13) 

(2.14) 



where we decomposed the partial derivative d'^ — u^^D -\- 
into "longitudinal" and "transverse" components 
D = u^di, and V = A^''^^^,, which in the local fluid 
rest frame reduce to the time derivative / = Df and 
spatial gradient V/. The first two equations describe 
the dilution of the local baryon and energy density due 
to the local expansion rate 0, while the third describes 
the acceleration of the fluid by the spatial (in the local 
frame) pressure gradients, with the enthalpy e-\-p act- 
ing as inertia. The 5 equations l|2.12|) - (|2.14|) for the 6 
unknown functions n, e, p, (remember that u^Ufi = 1) 
must be closed by supplying an Equation of State (EOS) 
p=p{e,n). 



III. NON-IDEAL FLUID DECOMPOSITION 

As the hydrodynamic evolution changes the local en- 
ergy and baryon density, microscopic processes attempt 
to readjust the local phase-space distribution to corre- 
sponding new local temperatures and chemical poten- 
tials. If this does not happen fast enough, the phase- 
space distribution will start to deviate from its local equi- 
librium form H2.4|l : f{x,p) = fcq(p-u{x);T{x), fi{x)) + 
Sf{x,p). The optimal values for the (readjusted) local 
temperature and chemical potential in the first term are 
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fixed by imposing the "Landau matching conditions" 'i'l 



E 



iu-prSfix,p)^0, (3.1) 



5m= (u-p) Sfix,p) - 0. (3.2) 

The remaining deviations Sf from local equilibrium gen- 
erate additional terms in the decompositions of N^, T^^ , 
and S*^: 



eq 



'5'cq 



(3.3) 

(3.4) 
(3.5) 



The new terms describe a baryon flow = A^'^N^, in the 
local rest frame, an energy flow = ^-^V^ + (where 
is the "heat flow vector") in the local rest frame, 
the viscous bulk pressure 11 = — iAp^T'*'' — p (which con- 
tributes to the trace of the energy momentum tensor), 
the traceless viscous shear pressure tensor tt^'^ — T^^'^^ = 
[ 1 (A'^'^A'^^-f A'^'^A^^) - iA^'^A'"^] Tra (where the ex- 
pression (/ij/) is a shorthand for "traceless and transverse 
to and defined by the projector in square 

brackets), and an entropy flow vector in the local 
rest frame. 

Inserting the decompositions 1)3. I|3.4|l into the con- 
servation laws (|2.1|l . (|2.2|l yields the non-ideal fluid equa- 
tions 



e = -{e+p+Il)e + TT^,^V 



(3.6) 

V-W + 2W-U, (3.7) 



(e+p+n) u^' = V^(p+n) - A^^Vvr^^ + tt^'^u 
Af'^'W^ + w^'e + {W-V)u 



(3.8) 



The matching conditions H3.1|l leave the choice of the 
local rest frame velocity it^ ambiguous. This ambigu- 
ity can be used to eliminate either V^^ from Eq. (|3.3(l 
("Eckart frame" , no baryon flow in the local rest frame 
Q), in which case the energy flow reduces to the heat flow 
vector or from Eq. ("Landau frame", 

no energy flow in the local rest frame ,4J), in which case 



there is a non-zero baryon flow ■ 



e+p 



q^ due to heat 



conduction in the local rest frame. (Intermediate frames 
are also possible, but yield no practical advantage.) For 
systems with vanishing net baryon number (as approxi- 
mately realized in RHIC collisions) the Eckart frame is 
ill-defined [31, so we will use the Landau frame. In this 
frame, for baryon-free systems with n = and no heat 
conduction, the non- ideal fluid equations H3.6|l - ()3.8|l sim- 
plify to 

e = -(e+p+n)6l + 7r^^V<'^M''\ (3.9) 
(e+p+n) = V^(p+n) - A^'^V^TT,^ -f ^''"it,. (3.10) 

The non-equilibrium decompositions H3.3|l - H3.5|l in- 
volve H-3-|-5=9 additional dynamical quantities, the 



"dissipative flows" 11, qf^, and tt'^" (the counting re- 
flects their transversality to u^^ and the tracelessness of 
tt'^'^). This means that we need 9 additional dynamical 
equations which should be compatible with the under- 
lying transport theory for the non-equilibrium deviation 
df{x,p). For the baryon-free case without heat conduc- 
tion, the number of needed additional equations reduces 
to 6. 



IV. KINETIC EQUATIONS FOR THE 
DISSIPATIVE FLOWS 



The key property of the kinetic equation governing 
the evolution of the phase-space distribution function 
/ = /cq+<5/ is that the collision term satisfies the sec- 
ond law of thermodynamics H2.3|l . i.e. entropy is pro- 
duced until the system has reached a new state of local 
thermal equilibrium. We don't want to solve the kinetic 
theory; instead, we want to write down a phenomeno- 
logical macroscopic theory which is consistent with the 
constraints arising from the underlying kinetic theory, in 
particular the 2"^^ law. The macroscopic theory will be 
constructed from an expansion of the entropy production 
rate in terms of the dissipative flows which themselves 
are proportional to the off-equilibrium deviation Sf of 
the phase-space distribution p. Assuming the latter to 
be small, |^/|^|/eq|, this expansion will be truncated 
at some low order in the dissipative flows SN^^, ST^^'^ . 
The expansion will involve phenomenological expansion 
coefficients which, in principle, should be matched to the 
kinetic theory (5|. In practice, they will often be con- 
sidered as phenomenological parameters to be adjusted 
to experimental data. In the end, the extracted values 
must then be checked for consistency with the entire ap- 
proach, by making sure that the dissipative corrections 
are indeed sufficiently small to justify truncation of the 
expansion a posteriori. 

The equilibrium identity H2.11|l can be rewritten as 



cq 



cq 



(4.1) 



where a=^, /3= ^, and l3i,= ^. The most general off- 
equilibrium generalization of this is 



(4.2) 

= p{a, I3)f3>' - aN^" + p^T^i" + Q^iSNt", ST""), 

where, in addition to the first order contributions im- 
plicit in the second and third terms of the r.h.s., in- 
cludes terms which are second and higher order in the 
dissipative flows 5N^ and 6T^^ . [Note that, by using 
the identity H2.11|l between the equilibrium quantities, 
Eq. H4.2I) can be written in the simpler-looking form 

the entropy production rate 

The form of the expansion (|4.2|l is constrained by the 
2"'^ law dfj.S^ >0. To evaluate this constraint it is useful 



^-\-Q'^ but this is not helpful for calculating 
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to rewrite the Gibbs-Duhem relation dp = s dT + ndfj, as 



B. Second order Israel-Stewart theory 



(4.3) 



With additional help from the conservation laws (|2. 1(1 and 
()2.2|l . the entropy production then becomes 



A causal theory of dissipative relativistic fluid dynam- 
ics is obtained by keeping up to terms which are sec- 
ond order in the irreversible flows. For simplicity we here 
consider only the baryon-free case n = = 0; for a gen- 
eral treatment see One writes 



(4.4) 



Using Eqs. (|3.3I3.4|I to express 5N^^ and ST^^" in terms of 
the scalar, vector and tensor dissipative flows 11, q'^, and 
TT^'', and introducing corresponding scalar, vector and 
tensor thermodynamic forces (in terms of gradients of the 
thermodynamic equilibrium variables) which drive these 
dissipative flows, X = —6 = —V -u^ X'" = ^-t^ — = 
-f^^ V''(^), andX^'^EEV^^u''^ (note that Xt'^'^X'^f"'^ 
is traceless and transverse to u), the 2"*^ law constraint 
can be further recast into 



(4.8) 



(with phenomenological expansion coefRcients /3o, P2) 
and computes (after some algebra using similar tech- 
niques as before) the entropy production rate as 



Td-S = n 



-9 - f3otl - TlTd^ 



2T 



(4.9) 



Tq/3 



2T 



Tdfj,S^ — UX q^X^ + TT^ X^y + Td^Q'^ > 0. (4.5) From the expressions in the square brackets we see that 



Note that the first three terms on the r.h.s. are first order 
while the last term is higher order in the dissipative flows. 



A. Standard dissipative fluid dynamics (flrst order 
theory) 

The standard approach (which can be found, for ex- 
ample, in 0) one neglects the higher order contributions 
and sets = 0. The inequality l|4.5|l can then always be 
satisfied by postulating linear relationships between the 
dissipative flows and the thermodynamic forces. 



the thermodynamic forces —9 and V 



are now self- 



consistently modified by terms involving the time deriva- 
tives (in the local rest frame) of the irreversible flows H, 
TTaf}- This leads to dynamical ("transport") equations for 
the latter. We can ensure the 2"^^ law of thermodynam- 
ics by again writing the entropy production rate in the 
form 1(4. 7|l (without the middle term), which amounts to 
postulating 



n 



1 



11 + (9 + ncr9. 



2(T 



--[U + C0], 



(4.10) 



n = -c^. 



e+p \T 



(4.6a) 
(4.6b) 
(4.6c) 



Taf3 



(4.11) 



with positive transport coefficients C > (bulk viscosity) , 
A > (heat conductivity) , and r/ > (shear viscosity) : 



Td-S 



n2 

C 



2AT 



2r, 



> 0. 



(4.7) 



(The minus sign in front of the second term is necessary 
because g^, being orthogonal to u^, is spacelike, q^ < 0.) 
These are the desired 9 equations for the dissipative flows. 

Unfortunately, using these relations in the hydrody- 
namic equations (|3.6|) - (|3.8|) leads to hydrodynamic evo- 
lution with acausal signal propagation: if in a given fluid 
cell at a certain time a thermodynamic force happens to 
vanish, the corresponding dissipative flow also stops in- 
stantaneously. This contradicts the fact that the flows re- 
sult from the forces through microscopic scattering which 
involves relaxation on a finite albeit short kinetic time 
scale. To avoid this type of acausal behaviour one must 
keep Q'^. 



Here we replaced the coefficients /?o,2 by the relaxation 
times = (^Po and = 277/32- In principle both C, 77 and 

, should be calculated from the underlying kinetic 
theory. We will use them as phenomenological param- 
eters, noting that for consistency the microscopic relax- 
ation rates should be much larger than the local hydro- 
dynamic expansion rate, ^9 ^1. 

Let us shortly discuss the approximation in the second 
equalities in Eqs. H4.10|l and (|4.11|) : We are using an ex- 
pansion scheme for the entropy production rate in which 
the thermodynamic forces and irreversible flows are as- 
sumed to be small perturbations. The approximation 
in Eqs. H4.1()|l . ((4.11|l neglects terms which are products 
of the irreversible flows with gradients of the thermody- 
namic equilibrium quantities which are of the same order 
as the thermodynamic forces. These terms are thus effec- 
tively of second order in small quantities and should, for 
consistency, be neglected relative to the other terms in 
the square brackets which are of first order. If one wants 



5 



to keep them (as done by Muronga 0,13), one should also 
keep third-order terms in the entropy flow vector Q'' for 
consistency. Of course, where the thermodynamic forces 
and irreversible flows are really small, it shouldn't matter 
whether we keep or drop these terms. In practice, how- 
ever, one will use this approach when dissipative effects 
are expected to be significant, and the dropped terms 
may not be extremely small. In this case we believe that 
dropping them is more consistent than keeping them. 

There is another reason for dropping these terms: 
without them, Eqs. (|4.1U|) . H4.11|l are relaxation equa- 
tions which describe (in the local rest frame) exponential 
relaxation (on the time scales „ ) of the irreversible 
flows to the values given by Eqs. (|4.6|1 in the first order 
theory. However, if these terms are kept, one has instead 
equations of the form 



n 

l+7nC 



CO + nC7n 
C 



(4.12) 



n 



l+7nC 



and similarly for the shear pressure tensor. One sees that 
both the kinetic relaxation time and the viscosity are 

^^j, 1 which involves 

the macroscopic expansion rate df^u^. This contradicts 
the intuitive expectation that these transport coefficients 
should be expressible through integrals of the kinetic col- 
lision term which involve only microscopic physics (cross 
sections, local densities, etc.) 

In the second order Israel-Stewart formalism, one 
thus solves the dissipative hydrodynamic equations 
(|3.8() simultaneously with the kinetic relaxation equa- 
tions (|4.1Q|I . H4.11|l for the irreversible flows. Let us now 
look at these equations in more detail when expressed in 
a global coordinate system (and not in local rest frame 
coordinates as done up to now). 



V. TRANSVERSE EXPANSION DYNAMICS IN 
SYSTEMS WITH LONGITUDINAL BOOST 
INVARIANCE 

We are restricting our discussion to systems with lon- 
gitudinal boost invariance. With this approximation we 
can describe the transverse expansion in very high en- 
ergy heavy-ion collisions in a domain near midrapid- 
ity. Boost-invariant systems are conveniently described 
in [T,x,y,rf) coordinates where r = \Jt'^—z'^ is longitu- 
dinal proper time, ?7 = ^ ln[(i+z)/(t— z)] is space-time 
rapidity, and r j\_ = (x, y) are the usual Cartesian coor- 
dinates in the plane transverse to the beam direction 
z. Boost-invariant systems are then characterized by 
macroscopic observables which are independent of i] and 
by phase-space distributions which depend only on the 
difference Y—ri (where F = i ln[(i?-|-p^)/(£'— p^)] is the 
momentum-space rapidity of a particle with longitudinal 
momentum and energy E) . We denote 2-dimensional 



vectors in the transverse plane by r^, dx, Vj_, etc., and 
generally use lowercase latin letters to denote vector and 
tensor components in this curvilinear space-time coordi- 
nate system. The metric tensor in this coordinate system 
reads 

= diag(l,-l,-l,-l/r2), 
5m„ = diag (1, -1, -1, -r^). (5.1) 

The flow velocity is parametrized as 

M™-7i(l,«„i;j„0) =7i(l,t;i,0) (5.2) 

where 

1 1 



7-L 



(5.3) 



with vanishing flow component in rj direction and 
transverse flow velocity v±{T,r±). For vectors and ten- 
sors, the usual Cartesian derivatives 9^ must be replaced 
by covariant derivatives, denoted by semicolons: 

= + (5.4a) 

d^T"^ ^ T"^,„ = a„r"' + n.^r'^' + T"'=rL, (5.4b) 

where V]^ = \g"^ {djgkrn+dkgmj"dm9jk) are the Chri- 
stoffel symbols. The only nonvanishing components of 
rj-fc are 



yv — yv 

J]T Til 



T 



(5.5) 



The time derivative in the local comoving frame and the 
local expansion rate are thus computed as 



D = u ■ d = '-f±[dT + v±-'V ±) 



(5.6) 



1 



e = d-u=-dr (r7_L) + Vi- (7^-^^) . (5.7) 

T 

If the expanding system has additionally azimuthal 
symmetry around the beam direction (for example, cen- 
tral collisions between spherically symmetric nuclei) , it is 
advantageous to replace the Cartesian transverse coordi- 
nates (x, y) by polar coordinates (r, 0) since macroscopic 
quantities are then (/)-independent: 

a;'" = (T,r,0,r?), (5.8a) 
5™" = diag (1, -1, -ly, -l/r^), (5.8b) 
5„„ = diag (1,-1, -r^-T^). (5.8c) 

This leads to the following additional (to Eq. (|5.5|l ') non- 
vanishing Christoffel symbols 



rf —rf — 4-1 r'' — —r 
The flow velocity now simplifies to 

1 

= 7r(l,i'r,0,0) with 7r 



(5.9) 



(5.10) 
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with radial transverse flow velocity v± =Vr{T, r) and 
vanishing flow components u*^ and u^. Correspondingly 
the time derivative in the local comoving frame and the 
local expansion rate reduce to 



D = U ■ d — "Ir [dr + Vrdr^ , 



= d ■ U = — dr (j^r 



1 



(5.11) 
(5.12) 



We will now treat the azimuthally symmetric and non- 
symmetric cases separately. 



A. (l+l)-dimensional viscous hydrodynamics: the 
azimuthally symmetric case 

Due to azimuthal symmetry and longitudinal boost in- 
variance, the n = (f> and n = r] components of the equa- 
tions of motion T™";™ =0 are redundant. Using the re- 
sults of Appendix^ the n — r and n~r components can 
be written as 

p + n + t'^ttI'' 



p + n + r^TT"^ 



(5.13) 



(5.14) 



With the shorthand notations T'"" 

T 
T 



and Vr — — these are brought into "standard 



rrT"", Vr = rTVr, 

- these are brouj 

(Cartesian) form" 



(5.15) 



Vr]=T 



(p + U 



The corresponding transport equations for the dissipative 
fluxes read (using the explicit expressions HA.11|I for the 



J 



shear tensor from Appendix Q 

1 



(dr+Vrdr'^TT''"'' 
dr + Wr^r In = 



1 



(e 7r 



.3 T 

2ri f 6 ^rVr 



(5.16) 



Similar equations were derived in |g (with extra terms, 
however, resulting from the higher order corrections in 
Eqs. (|4.1()|l and (|4.11|l which we argued should be ne- 
glected). These equations can be solved with the code 
LCPF CT [3, using subroutine LCPFCT for equ ation s 

and subroutine CNVFCT for equations 
First attempts at a numerical solution have been reported 
in Refs. 0, |^ , but a number of open questions remain at 
this point, and we leave a detailed discussion of the nu- 
merical aspects to a subsequent publication. 

The hydrodynamic equations require the equation of 
state p{e) for closure, i.e. after each transport step in 
time we must extract at each spatial grid point the boost 
velocity Vj. between the global and local rest frames and 
the local energy density e from the dynamical variables 
T"^"^ and T'^''. Equations IjA.Sp give the energy density as 



VrT^ 



(5.17) 



where the radial velocity Vr must be extracted from the 
implicit equation 



T^^ + p(e=T^^-w^r^'-) -f n - r^TT-^-^ - t^tt')'' 

(5.18) 

by a one-dimensional zero search. This is still the same 
degree of numerical complexity as in the ideal fluid case 
jlOj; we will see in the next subsection, however, that 
this part of the problem becomes numerically more in- 
volved for dissipative hydrodynamics without azimuthal 
symmetry. 



B. (2+l)-dimensional viscous hydrodynamics with longitudinal boost invariance 

In the absence of azimuthal symmetry, only the n = 77 component of the equations of motion r'"";m = is redundant 
(due to boost invariance). Using the relations HA.21|I from Appendix F^l the n — T,x,y components can be written as 



1/ \ D + n + T'^Tr'- 

-drirT^^j + 9^T^^ + dyry ^ - — — ^ , (5.19a) 

^^.(tT--) + 9,((r^^'-^^-)i.,) + dyi^ir^^TT^^vy') = -d^p+n+TT--) - dyTT-y, (5.i9b) 
^^.(tt-^) +d,(^{ry-TT-y)v^'^ +dy(^{ry-TT-y)vy^ = -d^n-y - dyip+n+nyy). (5.i9c) 

The further manipulation of these equations depends on our choice of independent shear pressure components as 
dynamical variables. In the following two subsections we explore two different choices, each with its own advantages 
and disadvantages. 
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1. Choosing n'^'^ , A = tt^^— tt^", and tt^'^ as independent dynamical variables 



In the first approach we select tt'^'^ , tt^^, and the difference A = 7r^^— tt^^ as independent dynamical components of 
the shear pressure tensor. The last choice has the advantage that it vanishes in the azimuthally symmetric case, thereby 
automatically reducing the number of independent dynamical variables. The choice of tt"^"^ instead of the orthogonal 
combination S = tt^^+tt^^ is a matter of taste and not essential. They are related by Eq. (IA.19p . tt'^'^ = S + t^tt''''. 

Introducing the shorthand notations y™" — -j-T™" (without the factor r this time), tt™" = ttt"'" , p — rp, II = Tn, 
Vi = — for i = x,y, and following the procedure of Appendix FSl equations Il5.19|l can be recast into "standard 
form" 

drf^-" + d^if-^'^v^) + dyif-^yvy) = -(p + n + t^tt'"'), (5.20a) 
drf-^ + a,((f-^-^^-)t;,) + dy(^{f^^~n-^)vy') = -d^ip+IL+n^n - dy^"" , (5.20b) 

drf^y + d^(^{f^y^n^y)v^^ + dy(^{f'-«-n''y)vy^ = -d^.^'^y - dyip+tv+f^yy), (5.20c) 



where we have resisted to insert the lengthy explicit ex- 
pressions (|A.22|) for tt"^"^, tt"^^, -n^y , tt^^, and tt^^. Note 
that the latter involve the velocities Vx and Vy, so in or- 
der to evaluate the sources at time step n which drive the 
propagation to time step n-l-1 we must explicitly solve for 
the velocities Vx^y{ri_) at time step n. We will return to 
this issue momentarily. 

The transport equations to be solved together with 
Eqs. (|5.20(l are (see Appendi J5|l 



[dr + Vxdx + Vydy^ vr'''' = 

(dr + WySy^Tr^^ : 

{dr + Vxdx + Vydy^ A = 

dr + Vrdr^H — 



1 

1 



u + ce 



277cr"" 

(5.21a) 
TT^^ - 2?7cr^ 

(5.21b) 
A - 2?7(t'^) , (5.21c) 

(5.21d) 



where a^^ , cr'^'^ , and a'^ are given by Eqs. l|A.23p . 

In order to compute the pressure p from the equa- 
tion of state p(e) , we calculate from the dynamical vari- 
ables T"^"^, T'^^, and T'^y the energy density by combining 
Eqs. (|A.21a|) - (|A.21cp : 



(Note that all viscous pressures cancel in this relation.) 
This requires the velocities Vx and Vy which are given 
(implicitly!) by 



(5.23a) 
(5.23b) 



T^^ +p(e) + n - TT^^ ■ 

T^^ -Kp(e) + n - TT^^' 



Since the vectors (T'^^^T'^y) and (tt"^^, tt'^^) are not par- 
allel to each other, the direction of the flow veloc- 
ity vj_ — (vx,Vy) is no longer given by the direction of 
(rprx ^rpry-^ as is the case in ideal fluid dynamics [loj . 
and Eqs. (|5.23l) can no longer be reduced to a simple 
one-dimensional zero search. Instead one must simulta- 
neously iterate two equations, one for the magnitude of 
the transverse velocity, 



^ - TT^^ +p{e) + n 



(5.24) 



VxT 



vyT^y. 



(5.22) 



(which is easily verified from Eq. (|A.21al) '). and one 
for the azimuthal angle of the velocity vector, (j)^ — 

tan~^ (^) ■ "^^^ latter is obtained by writing Vx — 

v±cos,(t)y = v±z, Vy — v± sin(/)i, = v±y^l—z'^, dividing 
the two equations (|5.23|l by each other and inserting 
Eqs. (|A.22a|) . IIA.22b|) : 



2v^x/l~z^T^y 



TT^^ _ !l J'(7r^^-r27r'"')(2z2-l) + ^ 



= 0. 



(5.25) 



Since Eq. H5.24|l requires knowledge of the direction of v±^ on the right hand side because of Eq. (|5.22|) . and 
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Eq. H5.25|l requires knowledge of v± , these two equations 
cannot be decoupled, and the iteration problem is gen- 
uinely 2-dimensional. This is of serious concern since this 
problem must be solved at every spatial grid point after 
each time step which makes it numerically very expen- 
sive. 



2. Choosing tt'^'^, tt^^, tt'^^, and tt^^ as independent 
dynamical variables 

This problem can be avoided if, instead of A, tt"^^ and 
n'^y are kept as dynamical variables which are directly 
evolved in time via their own kinetic transport equations. 
Defining the 2-dimensional vector M={Mx,My) = 
^rpTx_^Tx ^ rpTy_^Ty^ (thls Is just thc transverse momen- 
tum density vector without the shear pressure contribu- 
tions), we see from Eqs. (|5.23|l that v±^ and M are par- 
allel, va_-M = v^M with M = ^Ml+M^. Introducing 

further more A fp = T^^-tt^^, Eq. can with the help 

of Eq. (|A.20cp be rewritten as 

Ah-v^-M = M(i-v^M (5.26) 

which requires knowledge of only the magnitude of 
v±. The latter can then be obtained by a normal 1- 
dimensional zero search from Eq. (|5.24() which can also 
be written as 

- ^ 
'"^ ~ Mo+p{e=Mo-v^M) + W' ^^'^^^ 

and the velocity components are reconstructed from 

M 

Vx = VA^j^. (5.28) 



M 



Note that this procedure requires direct knowledge of 
^TT^ ^Tx^ g^^^ ^Ty g^-j- spatial grid points at each time 
step, so tt"^^ and tt'^^ cannot be computed from the con- 
straints IIA.22a|l . ljA.22bl) (for which the velocities and 
Vy would need to be already known). On the other 
hand, tt'^'^, tt"^^, and tt"^^ are not independent, but related 
by Eq. (|A.20cp . Since the suggested procedure requires 
propagating all three of these shear pressure components 
independently via kinetic transport equations (which, of 
course, should accurately preserve the constraint HA.20c|) 
if correctly implemented numerically), we have to solve 
one more kinetic transport equation (involving one phys- 
ically redundant component) than in the procedure of 
Section FVB II Instead of the three independent kinetic 
transport equations (|5.21a|l - H5.21c|) we have to solve four 
equations for tt"^"^, tt"^^, vr"^^, and vr''''. 

The set of equations to be solved simultaneously in 
this approach is given by Eqs. (|5.2U|) . H5.21b ..b.d). plus 
the following two equations: 



[pr + + Vydy)TT''y 



(5.29) 



The shear tensor components required here are given in 
Appendix [3 Eqs. (|A.25|I . 

Is this approach more economical than the two-dim- 
ensional zero search from the previous subsection? We 
believe so. If a one-dimensional zero search requires N 
iterations, each with K algebraic manipulations, a two- 
dimensional zero-search would require 0{N^-K-K') alge- 
braic manipulations at each spatial grid point and time 
step. Solving instead an additional kinetic transport 
equation for, say, tt"^^ requires 0{Ks-Kt) algebraic ma- 
nipulations at each time and grid point, where Ks is the 
number of algebraic steps required to evaluate the source 
cr'^^ and Kt is the number of algebraic steps involved in 
the time evolution algorithm. Taking the number K' of 
manipulations required to evaluate Eq. (|5.25|) to be com- 
parable to Ks^ solving the extra kinetic transport equa- 
tion should be numerically less expensive since we expect 
Kt to be significantly smaller than N^-K. The actual nu- 
merical implementation will tell whether this expectation 
is borne out. In any case, dissipative hydrodynamics is 
considerably more expensive than ideal fluid dynamics, 
and efficient coding will be required. 



VI. SUMMARY 

In this paper we derived explicit equations of motion, 
in a form that makes them directly amenable to publicly 
available transport algorithms for a causal theory 
of dissipative hydrodynamic evolution for relativistic vis- 
cous fluids such as those created in relativistic heavy-ion 
collisions. In doing so we followed the pioneering work 
by Israel and Stewart j5(| which was recently brought 
to wider attention and worked out in greater detail by 
Muronga |^. Our treatment is still not completely gen- 
eral in that it continues to assume boost-invariant ex- 
pansion along the beam direction (thereby reducing the 
spatial dimensionality of the problem by one) , but it goes 
beyond the existing literature 0, by allowing for ar- 
bitrary transverse expansion, without the additional re- 
striction of azimuthal symmetry around the beam direc- 
tion. It thus provides thc formal basis for a numerical 
calculation of elliptic flow in relativistic heavy-ion colli- 
sions including dissipative effects. Such calculations will 
be needed for the phenomenological determination of the 
viscosity of the quark-gluon plasma from heavy-ion col- 
lision data. 

We have also added to the discussion in 8] of az- 
imuthally symmetric (l-)-l)-dimensional viscous hydro- 
dynamics by reanalyzing the Israel-Stewart approach 
as presented in Ref. and pointing out a systematic is- 
sue with the expansion of the entropy current in terms of 
higher order terms in thc dissipative fluxes. The practical 
importance of the improved truncation scheme suggested 
here will only be fully assessed once the numerical imple- 
mentation of our equations Q has been thoroughly tested 
and becomes available for systematic investigations. 

The simplification of the problem resulting from az- 
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imuthal symmetry, by exploiting polar coordinates in the 
transverse plane, is significant. The effects of viscosity 
can be subsumed into an effective radial pressure, leav- 
ing the diagonal structure of the energy-momentum ten- 
sor in the local rest frame intact. As a consequence, the 
numerically critical problem of extracting at each time 
step from the dynamical components of T''" the local flow 
velocity and energy density in order to compute the pres- 
sure from the equation of state remains one-dimensional, 
i.e. of the same complexity as for ideal fluids. The addi- 
tional complexity resulting from dissipation thus resides 
entirely in the need for solving, together with the two hy- 
drodynamic evolution equations, three additional kinetic 
transport equations for the bulk viscous pressure and for 
two components of the shear viscous pressure. 

For the general situation without azimuthal symme- 
try, the use of polar coordinates (with their coordinate 
singularity at rj_ — 0) only complicates matters. It pro- 
vides no help towards solving the now in general two- 
dimensional selfconsistency problem associated with the 
extraction of the flow velocity and local energy density 
from the dynamical variables. We therefore use Carte- 
sian coordinates in the transverse plane, as has been the 
tradition in (2-|-l)-dimensional ideal fluid dynamics. Un- 
fortunately, this choice eliminates the possibility of cod- 
ing the equations in such a way that the code would au- 
tomatically take full advantage of all the simplifications 
resulting from azimuthal symmetry when handed an az- 
imuthally symmetric problem. Azimuthally symmet- 
ric (l-l-l)-dimensional expansion and asymmetric (2-1-1)- 
dimensional expansion require differently optimized algo- 
rithms. 

Using Cartesian transverse coordinates, we found a 
nice way of avoiding the above-mentioned two-dimen- 
sional nature of the numerically critical iteration prob- 
lem for the local energy density, by increasing the set of 
kinetic transport equations not by 1, but by 2 relative 
to the azimuthally symmetric case. By keeping one of 
the redundant components of the shear pressure tensor 
as a dynamical variable, we can again bring the iteration 
problem for the local energy density into scalar form. 
We believe that the expense for solving an additional 
transport equation, although not negligible, is less than 
that required for coping with a two-dimensional iteration 
problem at each time step at all spatial grid points. 

Compared to (2-f l)-dimensional ideal fluid dynamics, 
dissipative dynamics generates more complicated source 
terms for the three independent hydrodynamic evolution 
equations and requires the additional simultaneous solu- 
tion of five kinetic transport equations, one for the bulk 
viscous pressure and four for shear viscous pressure com- 
ponents (one of them being physically, but not algorith- 
mically redundant). Altogether the resulting increase in 
numerical complexity (compared to the ideal fluid case) 
is probably less than an order of magnitude. Given the 
increase in computer speed and power during the past 
decade, this should be manageable. 
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APPENDIX: SHEAR TENSOR AND VISCOUS 
PRESSURE TENSOR COMPONENTS 

1. Azimuthally symmetric systems 

For azimuthally symmetric systems we use polar co- 
ordinates in the transverse plane: a;™ = (t, r, 0, 7y) . In 
this coordinate system, the global frame and the local 
fluid rest frame are connected by a radial boost with 
velocity v = Vr (r, r) . Introducing the fluid rapidi- 
ty Ur = tanh~^ Vr such that 7^ = coshj/r, IrVr = sinhwr, 
the corresponding Lorentz transformation matrix is given 
by 



A"„(«,.) 



/ cosh i/r sinh j/,. 

sinh^r coshj/r 

10 

V 1 



(A.l) 



The projector transverse to the flow vector u"^ takes the 
form 



A" 



/-7>. 

V 









(A.2) 



!)■ 



(note that this differs from 

Due to azimuthal symmetry and longitudinal boost- 
invariance, all mixed components involving indices or 
?; of the stress tensor 



cr"" EE V<"u"> 

and the shear pressure tensor tt™" vanish: 



j-4>V 



r'lr 



0, 
0. 



(A.3) 



(A.4) 



This leaves only two independent components for the 
shear pressure tensor which is constrained by the con- 
ditions of tracelessness 



T TT ' ' 



and of orthogonality to =7^(1, Wr, 0, 0): 



VrT^ 



(A.5) 



(A.6) 



Equations (|A.6|I result from the n = T and r components 
of UmT^™"^ — 0; the other two components yield redundant 
equations. Equations (|A.5p . (|A.6p can be combined to 
yield 



-7^(r"7r^*-hr2^''") 



(A.7) 
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Identical relations hold between the corresponding com- 
ponents of the shear tensor ct™". Equations HA.5() and 
(|A.7|I can be used to eliminate tt"^"^, tt'"^ and n^^ in favor 
of tt"^-^ and tt't^ which we keep as independent dynamical 
variables. Equation IjA.ll) shows that these components 
are not affected by the radial boost, i.e. the (p(j) and rji] 
components of the shear pressure tensor are identical in 
the global frame and in the local fluid rest frame. The 
same is not true for the other non-zero components of 
the shear pressure tensor. 

This observation leads to a particularly simple struc- 
ture of the energy-momentum tensor in the global frame, 
including viscous pressure terms. Using Eqs. (|A.5p - HA.7p 
in the expressions for T'^'^ and T'^^ (which are the two 
independent components being evolved by the hydrody- 
namic evolution equations) we find 



2. Systems without azimuthal symmetry 

Without azimuthal symmetry, there is no advantage 
in using transverse polar coordinates, and it is simpler 
to use Cartesian coordinates in the transverse plane: 
= (r, X, y, if). The global frame and the local fluid rest 
frame are now related by a Lorentz boost with velocity 
v±=Vxix,y,T) + Vy{x,y,T) By = {vx,Vy), described 
by the Lorentz transformation matrix 



(e+p+n)72 



(p+n) 



(A.8a) 



/ 7± ll-Vx ^ l±Vy \ 

l±v., l + (7^-l)^ (7±-l)^ 

l±vy h±-l)^ l + (7±-l)f 
\ ^ ^ l/ 

(A.12) 

The projector transverse to the flow vector u™ takes the 
form 



T^"- = {e+p+n)-ff.Vr+TT^'' ^ {e + VrhrVr, (A.8b) 










VxVy 


where 


A 


m _ 

n ~ 




~Vy 















Vr ^p + n-r^n''"'' -t^tt'^'K (A.9) 










Equations (jA.Sp are the same expressions as for the ideal 


In the local 


rest 


frame 


the e 


reads 










fluid, except for the replacement of p by the effective 


( e 








radial pressure Vr which has no explicit dependence on 










' \ 


the flow velocity Vr- 


rpmn 













In the local rest frame, the energy-momentum tensor 








P+n 





is diagonal: 















Vy 

VxVy 









(A.13) 











\ 





Vr 














p+n+r^TT'l"!' 





Vo 








/ 



(A.IO) 



The non-vanishing components of the shear tensor 
^mn _y{m^n) g^^.^ found (aftcr some algebra which prop- 
erly takes into account the Christoffel symbol contribu- 
tions to the covariant derivatives) to be 



/O 

^xx ^xy Q 
^xy ^yy q 

V 71-'''' 

(A. 14) 

The shear pressure tensor tt'"'* is no longer diagonal. It 
may be useful to have expressions for the shear pressure 
tensor components tt"^" in the local rest frame in terms 
of the hydrodynamic solution for the energy-momentum 
tensor T'"" in the global frame. To this end we follow 
Ref. 8] and introduce the mutually orthogonal rest frame 
4- vectors 



1 



7r 

r 

^r^r 



-Ir 



20 



7r 

r 



IrVr 



7r o - drlr " 9r(7r«r) 



.3 

(1-7.') 



/ 20 ^ > _ 7r 

V 3 r 7 



VrU 



(A.lla) 
(A.Ub) 

(A.llc) 

(A.Ud) 
(A.Ue) 



with u-u = 

n 



u = 

r = 

1, i-i 



(1,0,0,0), 
(0,0,1,0), 



h" 



(0,1,0,0), 
: (0,0,0,1/t) 



(A.15) 



-J-J - 



-h-h- 



P 



TT 



T TT ' ' 



1 

3 

2 

3 
2 

3 
2 



(A.16a) 
(A.16b) 



Note that IjA.llbl) differs from the corresponding result 
(44) in 8] by a metric factor , and that the last expres- 
sion in (jA.llcp has an extra factor 7^ when compared to 
Eq. (43) in 8]. 



1, such that 

(i-T-t) + U-T-j) + (h-T-h) , 

p-U 

{i-T-i)~^-({rT-]) + {h-T-h)), (A.16c) 

(j-T-j)-^({z-T-i) + {h-T-h)), (A.16d) 

{h-T-h)~^-({i-T-i) + {rT-j)). (A.16e) 



The right hand sides are Lorentz invariant expressions 
and thus can be evaluated in any reference frame. The 
vectors u'", z*", j*", /i"' in the global frame are obtained 
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by applying the Lorentz boost ljA.12p to the rest frame The constraints from tracelessness, 
vectors (IA.15|I : 



(A.17a) 



TT =7r + TT"" + T TT ' ' 



J 



(7x^;.,l+(7±-l)%,(7±-l)^,0), (A.17b) 

2 

7±i;^,(7±-l)^,l+(7±-l)^,0), (A.17c) 

^1 ^ I 



and orthogonality to 



- 0,0,0, 



1 



(A.lTd) 



(A.19) 



(A.20a) 
(A.20b) 
(A.20c) 



From (|A.17d|) it follows immediately that the 7777 com- 
ponent of T"™ is identical in the global and fluid rest then leave us with three independent components of tt™". 
frames, i.e. that tt'''' = tt'''' . Longitudinal boost invari- Which should we select? From the above it is obvious 
ance implies that that ir^^ should be one of them. For the other two there 

are two different possibilities, following from two different 

^jjT _ ^rfx ^ ^r)ij ^ ^jjr _ ^7^x ^ ^771/ _ q_ chains of reasoning, as explained in Sections IV B II and 

(A.18) rVB2l 



Before giving the explicit expressions for the shear tensor needed in each case, let us first generalize the relations 
(jA.8|l . which are needed to bring the hydrodynamic equations into standard form: 



rpTT 






hn)7i - 


(p+n) + 








J-iTX 


= {e- 




^^hlvx 




[T^^+P- 






rpry 


= ie- 






+ n^y = 




^n-TT^^)Vy + 


TT-'y, 


rpXX 




^P- 




+ P + II 


j^^xx ^ 


[T--^^rx^V^ 


+ p + n + 7r^^, 




= ie- 


hpH 


y-iihlvx 


Vy + 




-Tr^'')Vy + TT^J' 


= {ry-T:^y)vx + TT^^ 




= {e- 






+ P + U 


+ Txyy = 


{T^y-T:-^y)vy 


+ p + ii + Tvyy , 



(A.21a) 
(A.21b) 
(A.21c) 
(A.21d) 
(A.21e) 
(A.21f) 



If we choose tt"^"^, tt'''' and A = 7r^^— tt"" as independent dynamical variables, we can solve the constraints (|A.19p . 
(jA.20|l to gives the following expressions for the dependent components of tt™": 



2z;^7r^^ = TT^M 1 



2 2 



2v^T:^y = 1 - 



(A.22a) 
(A.22b) 



2 2 



(A.22c) 

(A.22d) 
(A.22e) 



In this approach we need the following components of the stress tensor cr™" as source terms for the kinetic transport 
equations for vr'^^, vr'''', and A: 



1 



-(7i-i) + a,7±-x^(7i) 



tM3 



7± 



„xx „yy 
O = (J — O 



1(2 + il{vl-vl^ + dyi^^vy) - dxi^^vx) - Id{jI{vI-vI)) . 



(A.23a) 
(A.23b) 
(A.23c) 
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If instead of A we keep tt"^^, tt"^^ as dynamical variables (as explained in Section IV B 2\i . we don't need Eqs. (jA.22a|) 
and (|A.22b|l . but we still must express tt^^, tt^^, and tt^^ (which appear as sources on the right hand sides of 
Eqs. (|5.20bp and H5.2()c(l ^ in terms of those components for which we solve kinetic transport equations. To this end 
we rewrite the constraints (IA.19II . (|A.20p as 

TT^^ = l.(^y-^(n^^^r\'"^) + v^TT^^ - Vyn'-yy (A.24a) 

TT^y = -J- (v2(^rr_^2^„„) _^ ^^^^ry _ „^^r:r^ ^ (A.24b) 
TT^^ = --^(y^VyiTT^^-T^TT'"^) - W^Tt"^ - W^^TT^"^) . (A.24c) 

The shear tensor components required on the right hand sides of Eqs. H5.29|l are 

11 1 6 

^'"^ - -^a,7± + \dr{l^vy) - \D{^lvy) + ^ 7i«y (A.25b) 
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